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Abstract 

Background: Mosquito control programmes using chemical insecticides are increasingly threatened by the 
development of resistance. Such resistance can be the consequence of changes in proteins targeted by insecticides 
(target site mediated resistance), increased insecticide biodegradation (metabolic resistance), altered transport, 
sequestration or other mechanisms. As opposed to target site resistance, other mechanisms are far from being fully 
understood. Indeed, insecticide selection often affects a large number of genes and various biological processes 
can hypothetically confer resistance. In this context, the aim of the present study was to use RNA sequencing 
(RNA-seq) for comparing transcription level and polymorphism variations associated with adaptation to chemical 
insecticides in the mosquito Aedes aegypti. Biological materials consisted of a parental susceptible strain together 
with three child strains selected across multiple generations with three insecticides from different classes: the 
pyrethroid permethrin, the neonicotinoid imidacloprid and the carbamate propoxur. 

Results: After ten generations, insecticide-selected strains showed elevated resistance levels to the insecticides used 
for selection. RNA-seq data allowed detecting over 13,000 transcripts, of which 413 were differentially transcribed 
in insecticide-selected strains as compared to the susceptible strain. Among them, a significant enrichment of 
transcripts encoding cuticle proteins, transporters and enzymes was observed. Polymorphism analysis revealed over 
2500 SNPs showing > 50% allele frequency variations in insecticide-selected strains as compared to the susceptible 
strain, affecting over 1000 transcripts. Comparing gene transcription and polymorphism patterns revealed marked 
differences among strains. While imidacloprid selection was linked to the over transcription of many genes, 
permethrin selection was rather linked to polymorphism variations. Focusing on detoxification enzymes revealed 
that permethrin selection strongly affected the polymorphism of several transcripts encoding cytochrome P450 
monooxygenases likely involved in insecticide biodegradation. 

Conclusions: The present study confirmed the power of RNA-seq for identifying concomitantly quantitative and 
qualitative transcriptome changes associated with insecticide resistance in mosquitoes. Our results suggest that 
transcriptome modifications can be selected rapidly by insecticides and affect multiple biological functions. 
Previously neglected by molecular screenings, polymorphism variations of detoxification enzymes may play an 
important role in the adaptive response of mosquitoes to insecticides. 
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Background 

Mosquitoes are vectors of several human diseases repre- 
senting a major burden for public health worldwide [1]. 
Half of the world's population is exposed to malaria while 
dengue fever represents a burden in over 100 countries 
with 2.5 billion people at risk [2,3]- Since the 1950s, chem- 
ical insecticides have been massively used for controlling 
mosquito populations but their efficacy is now threatened 
by resistance mechanisms developed by insects. In absence 
of efficient alternatives, characterizing molecular mecha- 
nisms underlying resistance is a key step for improving re- 
sistance management strategies. 

Resistance to insecticides can be the consequence of dif- 
ferent mechanisms, such as a mutation of the proteins tar- 
geted by the insecticide (target-site insensitivity), a lower 
penetration of the insecticide, its sequestration, or its 
biodegradation (metabolic resistance) [4,5]. Target-site in- 
sensitivity and metabolic resistance are known as the two 
main resistance mechanisms in mosquitoes [4,6]. Muta- 
tions causing target-site insensitivity are well-characterized 
in mosquitoes and molecular tests detecting these muta- 
tions are available for species of public health importance 
[7-12]. Metabolic resistance has been reported worldwide 
and usually involves detoxification enzymes such as cyto- 
chrome P450 monooxygenases (P450s or CYPs for genes), 
carboxy/cholinesterases (CCEs), glutathione S -transferases 
(GSTs) and UDP glucosyl-transferases (UGTs) [4,6,13,14]. 
However, due to the large number of mosquito genes en- 
coding detoxification enzymes [15-17] pinpointing those 
responsible for resistance remains challenging [18]. 
Metabolic resistance has been mostly associated with an 
increased level of detoxification enzymes in resistant 
populations and multiple candidate genes have been 
identified by microarray screenings [6,18-21]. In con- 
trast, polymorphism variations potentially affecting the 
functionality of detoxification enzymes have been hardly 
studied in mosquitoes despite evidences suggesting that 
this phenomenon may play a role in insecticide resistance 
[14,22]. Recently, polymorphism of a P450 gene has been 
associated with pyrethroid resistance in the mosquito 
culex pipiens [23] and a reduction of sequence diversity in 
two P450 genes conferring resistance to pyrethroids has 



been observed in Anopheles funestus [24]. This suggests 
that a deep analysis of the polymorphism associated with 
resistance can improve our understanding of mechanisms 
developed by mosquitoes to resist insecticides. Today, 
such knowledge gap can be overcome by high throughput 
sequencing approaches such as mRNA sequencing (RNA- 
seq), which can generates concomitantly gene expression 
and polymorphism data over the whole transcriptome 
from a single experiment [25-27]. 

In this context, the aim of the present study was to use 
RNA-seq for investigating transcription level and poly- 
morphism variations associated with adaptation to three 
insecticides from distinct chemical families in the mos- 
quito Aedes aegypti. A susceptible strain was selected with 
the pyrethroid permethrin, the neonicotinoid imidacloprid 
or the carbamate propoxur, to produce three resistant 
strains. After ten generations of selection, the constitutive 
resistance level of each resistant strain was measured and 
the transcriptome of each strain was deep sequenced. 
After mapping cDNA reads to the genome, gene expres- 
sion and polymorphism variations linked to insecticide 
selection were identified and compared across strains. Re- 
sults are discussed in regards of known and new putative 
adaptive mechanisms conferring insecticide resistance in 
mosquitoes. 

Results 

Insecticide resistance levels 

After ten generations of larval selection with the insecti- 
cides permethrin (Perm-R strain), imidacloprid (Imida-R 
strain) or propoxur (Propo-R strain), bioassays revealed a 
constitutive increased resistance of each selected strain to 
its respective insecticide as compared to the parental sus- 
ceptible strain (Table 1). Resistance levels of the Perm-R 
and Imida-R strains were moderate but significant (3.78 
fold and 5.85 fold respectively). Although significant, the 
resistance level of the Propo-R strain to propoxur was 
considerably lower (1.84 fold). 

Sequencing, read mapping and genome re-annotation 

More than 269 million 75 bp cDNA reads were sequenced 
across all samples (Additional file 1: Table SI). Each 



Table 1 Resistance levels after insecticide selection 



Insecticide 


Strain 


LCso (ng/L) 


Cl95% (MQ/L) 


RR50 


Ci95% 




Susceptible 


2.52 


2.29 - 2.77 






Permethrin 












Perm-R 


947 


6.21 - 14.43 


3.78 


2.24 - 630 




Susceptible 


240.30 


208.74 - 276.67 






midacloprid 












Imida-R 


1406.40 


1236.24- 1599.87 


5.85 


4.47 - 7.66 




Susceptible 


441 .20 


392.16 -496.34 






Propoxur 










Propo-R 


813.60 


743.47 - 890.23 


1.84 


1.50 - 2.27 



Resistance ratios were computed from LC50 values as compared to the susceptible parental strain. 
Resistance ratios with confidence intervals not overlapping the value of 1 are shown in bold. 
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mosquito strain was represented by two cDNA libraries 
with an average of 33.6 million reads per library. More 
than 80% of reads were successfully mapped on the Ae. 
aegypti genome. Filtering on sequence quality and map- 
ping score retained more than 174 million reads (~ 65%). 
For each strain, RPKMs (Reads Per Kilobase of exon 
model per Million reads) obtained from the two cDNA li- 
brary replicates were well correlated (Additional file 2: 
Figure SI), indicating moderate variations between tech- 
nical replicates. In consequence, reads from technical rep- 
licates were pooled for further analyses. 

A transcription signal was detected for 85% of known 
Ae. aegypti genes. Comparing transcript coverage between 
strains revealed similar distributions with RPKMs span- 
ning more than 6 Logs and a median transcript coverage 
of 4.4 RPKM (-300 reads/Kb) (Additional file 3: Figure 
S2). Comparing genome annotation with the distribution 
and coverage of reads suggested incorrect gene boundary 
annotations for more than 3000 transcripts and identified 
more than 500 isolated novel transcription events (NTEs) 
based on their transcription signal, structure and high dis- 
tance to known transcripts (Additional file 4: Table S2 and 
Additional file 5: Table S3). Distribution of mapped 
reads across the whole Ae. aegypti genome can be 
accessed at http://vectorbase.org using the 'configure this 
page and 'RNAseq alignments' options of the genome 
browser (tracks 'Bora-Bora control! 'Perm-RJ 'Imida-R' and 
'Propo-R'). 

Differential transcription in insecticide-selected strains 

Differential transcription analysis was performed on the 
13105 transcripts showing a transcription signal higher 
than 0.5 RPKM in all strains. A total of 463 transcripts (~ 
3.5%) including 413 known transcripts and 50 NTEs were 
considered differentially transcribed in any insecticide- 
selected strain as compared to the susceptible strain (> 3 
fold in either direction and adjusted P value < 10"^^; 
Table 2). Such threshold appeared biologically relevant as 
less than 2% of transcription ratios belonging to 140 
housekeeping genes (all ribosomal proteins, actins and 



tubulins) were found differentially transcribed as com- 
pared to the parental strain. Cross-comparison of tran- 
scription ratios (TRs) obtained from RNA-seq and DNA 
microarray was performed for the Imida-R strain from the 
same biological samples (Additional file 6: Figure S3). This 
comparison revealed a good correlation between the two 
techniques (r = 0.86, slope ~ 1, most variations < 2 fold). 

The balance between over- and under transcription was 
contrasted between each strain (Figure 1). The Imida-R 
strain showed the widest response to insecticide selec- 
tion with 227 transcripts specifically over transcribed as 
compared to the susceptible strain. In contrast, fewer 
transcripts were affected in Perm-R and Propo-R strains 
with the majority of them being under transcribed. A total 
of 96 transcripts were found differentially expressed in 
multiple strains including 17 and 20 transcripts over- and 
under transcribed in all strains respectively. 

A clustering analysis based on TRs as compared to the 
susceptible strain was performed on the 413 known tran- 
scripts differentially transcribed in any insecticide-selected 
strain (Figure 2 and Additional file 7: Table S4). This ana- 
lysis confirmed the specific over transcription of several 
genes in the Imida-R strain and identified nine main tran- 
script clusters based on their expression profile across 
strains. Assigning known transcripts to biological categor- 
ies and comparing their frequency to all detected tran- 
scripts revealed protein families or biological functions 
enriched in the different clusters (Figure 2). When consid- 
ering all clusters as a whole, a significant enrichment of 
transcripts related to immunity and cuticular proteins was 
observed. Cluster 1, representing transcripts strongly over 
transcribed in Perm-R strain included 5 hexamerins asso- 
ciated to cellular trafficking. These hexamerins showed 
very high TRs in the Perm-R strain (from 10 to 150 fold) 
but were also over transcribed in other insecticide- 
selected strains (up to 12 fold). Cluster 2, representing 
transcripts over transcribed in all insecticide-selected 
strains, showed a significant enrichment in detoxification 
enzymes. Cluster 3, representing transcripts strongly under 
transcribed in the Perm-R strain but not in others, was 



Table 2 Differential transcription analysis overview 





Perm-R 




Imida-R 




Propo-R 




Any strains 






Transcripts 


% 


Transcripts 


% 


Transcripts 


% 


Transcripts 


% 


AC test P value 


2858 


21.8 


2942 


22.5 


1881 


144 


4188 


32.0 


AC test P value and FC >3 


181 


1,4 


339 


2.6 


131 


1.0 


463 


3.5 


Over transcribed 


34 


0.3 


259 


2.0 


38 


0.3 


279 


2.1 


Known transcripts 


34 


0.3 


235 


1.8 


37 


0.3 


255 


1.9 


New putative transcripts 


0 


0.0 


24 


0.2 


1 


0.0 


24 


0.2 


Under transcribed 


147 


1.1 


80 


0.6 


93 


0.7 


239 


1.8 


Known transcripts 


131 


1.0 


70 


0.5 


76 


0.6 


210 


1.6 


New putative transcripts 


16 


0.1 


10 


01 


17 


01 


29 


0.2 
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Over-transcribed versus susceptible strain Under-transcribed versus susceptible strain 

Figure 1 Transcripts differentially expressed after insecticide selection. For each Venn diagram section, tlie numbers of transcripts 
differentially expressed in any insecticide-selected strain as compared to the susceptible strain are indicated. 



significantly enriched in cuticle proteins. Clusters 4 and 5 
were characterized by transcripts showing a strong over 
transcription in the Imida-R strain and were significantiy 
enriched in cuticular proteins. A significant enrichment in 
proteins potentially involved in immunity was also ob- 
served for cluster 4. 

Polymorphism variations in insecticide-selected strains 

A total of 220,499 SNP loci were identified between any 
strain and the Ae. aegypti reference genome. As expected 
in populations under selection and experiencing recurrent 
bottlenecks, a notable lower variability in insecticide- 
selected strains was observed as compared to the parental 
susceptible strain (Table 3). Comparative polymorphism 
analysis identified several alleles differentially represented 
between insecticide-selected strains and the parental sus- 
ceptible strain (> 50% allele frequency difference, referred 
to as 'differential SNPs')- The Perm-R strain showed more 
differential SNPs (1315) as compared to other insecticide- 
selected strains (811 and 812 for the Imida-R and Propo-R 
strains, respectively). Predicted genie effects of these varia- 
tions were equally distributed among insecticide-selected 
strains with ~ 2.5% of them located in 5'UTRs, ~ 9% in 3' 
UTRs, ~ 60% within transcript coding sequences, ~ 5% in 
introns and ~ 20% within 1 kb of gene boundaries (re- 
ferred to as 'intergenic'). No significant correlation was 
observed between transcripts differentially transcribed 
and the presence of differential SNPs in their 5' or 3' 
UTRs (not shown). More than 1080 transcripts were af- 
fected by differential SNPs in at least one insecticide- 
selected strain (Figure 3). Up to 22 differential SNPs per 
transcript were observed, with 21 transcripts affected 
by > 10 differential SNPs. The Perm-R strain showed a 
higher number of transcripts affected by differential 
SNPs as compared to other strains (582 versus 415 and 
420). Only 25 transcripts were affected by differential 
SNPs in all insecticide-selected strains. No differential 
SNP was found in transcripts encoding known insecticide 



target proteins such as the para-type voltage-gated sodium 
channel (permethrin), the acetylcholinesterase (propoxur) 
or the nicotinic acetylcholine receptor (imidacloprid). 

More than 1650 differential SNPs were detected within 
coding regions (Figure 4 and Additional file 8: Table S5). 
Clustering analysis based on allele frequency variations 
between each insecticide-selected strain and the parental 
strain evidenced the differential response of each strain 
to insecticide selection. The Perm-R strain showed larger 
allele frequency variations as compared to other strains 
(median frequency variations of 50.2%; 29.8%; 27.9% for 
Perm-R, Imida-R and Propo-R respectively). When com- 
paring all detected transcripts with those affected by differ- 
ential SNPs (Figure 4, all clusters), a significant enrichment 
of transcripts encoding detoxification enzymes was de- 
tected. Such enrichment was also observed when consider- 
ing variations predicted as non- synonymous only. When 
considering each cluster independently, clusters 1 and 9, 
representing differential SNPs specific to the Imida-R 
strain, did not show any significant enrichment in any bio- 
logical category. Clusters 2, 4 and 6 were composed of dif- 
ferential SNPs specific to the Perm-R strain and revealed a 
significant enrichment in transcripts encoding detoxifi- 
cation enzymes, with multiple cytochrome P450s af- 
fected. As observed for transcription level variations, 
several transcripts encoding hexamerins were affected 
by large polymorphism variations in the Perm-R strain 
(Additional file 8: Table S5). Clusters 3 and 8, represent- 
ing differential SNPs found in the Imida-R and Propo-R 
strains, also showed an enrichment in detoxification en- 
zymes, as well as an enrichment in transcripts involved 
in immune response. 

Gene ontology terms enrichment analyses 

Gene Ontology (GO) term enrichment analyses were per- 
formed independently on transcripts significandy over- or 
under transcribed or affected by differential SNPs in their 
coding region (Additional file 9: Figure S4). These analyses 
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Transcription ratio (fold) 
+4 



r 



All detected 
transcripts 



Imida-R Perm-R Propo-R 





12473 
(6181) 



16 
(9) 



23 
(16) 



27 
(25) 



^1 



114 
(45) 



101 
(40) 



(8) 



79 
(34) 



(2) 



^ (23) 

413 

(202) 



■ Detoxification 

■ Kinases/Phosphatases 
Other enzymes All 

■ Cuticle clusters 

■ Immunity 
Hormones and neurotransmitters signalling 

■ Transcription factors 

_ Intra or extracellular trafficking/Chaperonin 

■ Structure 

Figure 2 Clustering of transcripts differentially expressed 
across strains. The analysis was performed on the 41 3 
known transcripts significantly differentially expressed in any 
insecticide-selected strain compared to the susceptible strain. 
Clustering was based on Euclidean distance of fold changes as 
compared to the susceptible strain and complete linkage algorithm. 
Pie charts describe biological functions affected within main clusters 
based on the number of transcripts assigned to each function. Stars 
indicate biological functions significantly enriched compared to their 
representation among all detected transcripts (Fisher's test adjusted 
P value: * < 0.05, ** < 0.01 , *** < 0.001). The total number of 
transcripts constituting each cluster is indicated. The number of 
transcripts with predicted functions used for building each pie chart 
is shown within brackets. 



Table 3 Differential polymorphism analysis overview 





Perm-R 


Imida-R 


Propo-R 




SNPs 


% 


SNPs 


% 


SNPs 


% 


Detected SNPs 


143737 


100 


143645 


100 


1 33892 


100 


Differential SNPs 


1315 


0.9 


811 


0.6 


812 


0.6 


Genie consequences* 


1403 


100 


881 


100 


908 


100 


5' UTR 


32 


2.3 


28 


32 


25 


2.8 


Synonymous coding 


715 


51.0 


442 


50.2 


453 


49.9 


Non-synonymous Coding 


137 


9.8 


95 


10.8 


108 


11.9 


Intronic 


66 


4.7 


45 


5.1 


54 


5.9 


3' UTR 


121 


8.6 


88 


10.0 


87 


9.6 


Intergenic 


325 


23.2 


184 


20.9 


181 


19.9 



^Consequences of allelic changes are only indicated for differential SNPs 
(i.e. alleles showing > 50% allelic frequency difference in any selected strain as 
compared to the susceptible strain}. These consequences were computed 
based on genome annotation. 



confirmed the distinct response of each strain to insecti- 
cide selection. Among GO terms representing genes over 
transcribed in insecticide-selected strains, terms related to 
cuticle proteins were strongly enriched in the Imida-R 
strain. Terms related to oxygen transport or to hexamer- 
ins were enriched in both the Perm-R and Imida-R strains, 
together with those related to iron homeostasis and iron 
binding. Among GO terms represented by under tran- 
scribed genes, multiple terms related to immune response 
were enriched in insecticide-selected strains. When con- 
sidering transcripts affected by differential SNPs, no GO 
term was enriched in the Imida-R strain while only 




Figure 3 Transcripts affected by differential SNPs. For each Venn 
diagram section, the number of transcripts affected by differential 
SNPs is shown. The number of transcripts affected by differential 
SNPs predicted as non-synonymous according to genome 
annotation are shown within brackets. 



David ef al. BMC Genomics 2014, 15:174 
http://www.biomedcentral.com/1471-2164/15/174 



Page 6 of 1 5 




Hormones and neurotransmitters signalling 
■ Transcription factors 

? Intra or extracellular trafficking/Chaperonin 
r Structure 

Figure 4 Clustering of differential SNPs across strains. The 

analysis was performed on all differential SNPs falling within coding 
regions. Clustering was based on Euclidean distance of allele 
frequency variations as compared to the susceptible strain and 
complete linkage algorithm. Pie charts describe biological functions 
affected within main clusters based on the number of differentia 
SNPs affecting transcripts assigned to each function. Stars indicate 
biological functions significantly enriched compared to their 
representation among all detected transcripts (Fisher's test adjusted 
P value: * < 0.05, ** < 0.01 , *** < 0.001). For each cluster, the total 
number of differential SNPs is indicated as plain text while those 
affecting transcripts of known function (used for pie charts) are 
shown within brackets. 



'dehydrogenase activity' was enriched in the Propo-R 
strain. In contrast, the Perm-R strain revealed an enrich- 
ment of all terms related to cytochrome P450s such as 
'monooxygenases activity', 'electron carrier activity', 'tetra- 
pyrrole binding', 'heme binding' and 'iron ion binding'. 



Terms related to hexamerins were also enriched in Perm- 
R together with terms related to various enzyme families. 

Focus on transcripts potentially involved in insecticide 
detoxification 

Focusing on cytochrome P450 monooxygenases (CYP 
genes. Figure 5) revealed that few CYPs showed signifi- 
cant transcription level variations, with none being over 
transcribed in the Perm-R strain and only four CYPs 
{CYP6BB2, CYP9M9, CYP6N9 and CYP6Z8) being over 
transcribed in the Imida-R and Propo-R strains. Four 
CYPs were under transcribed in the Perm-R strain 
{CYP6M5, CYP6M10, CYP6M9 and a CYP12F). This 
CYP12F was also under transcribed in the two other 
strains. When considering differential SNPs, the Perm-R 
strain carried a much higher number of differential 
SNPs affecting P450s as compared to other strains. The 
19 CYP genes specifically affected by permethrin selec- 
tion included 13 CYP6s belonging to a dense P450 
cluster located on supercontig 1.371 composed of six 
CYP6NS, five CYP6Ms, four CYP6Zs, CYP6S3 and 
CYP6Y3 (Additional file 10: Figure S5). In contrast, only 
5 and 4 CYPs were affected by differential SNPs in the 
Imida-R and Propo-R strains respectively. Several CYPs 
were affected by variations predicted as non-synonymous 
(26 variations affecting 16 genes) while 17 variations af- 
fecting 14 genes where located within substrate recogni- 
tion site regions (SRS), potentially affecting the active site 
of these P450s. 

When considering other transcripts potentially involved 
in resistance (Figure 6), the strong response of the Imida- 
R strain to insecticide selection through transcription level 
modifications was confirmed, with several detoxification 
enzymes being over expressed including the glutathione 
S-transferase GSTD4, one glycosyltransferase together 
with multiple oxidases, peroxidases and dehydrogenases. 
Several kinases were also specifically over-transcribed in 
the Imida-R strain. Fewer transcripts were affected in the 
Perm-R and Propo-R strains, including one aldo-keto re- 
ductase, one glycosyltransferase, one aldehyde oxidase and 
one alcohol dehydrogenase. As opposed to CYP genes, 
differential SNPs affecting other detoxification genes 
were well-balanced between strains with 66 variations 
affecting 26 distinct genes. Among them, ten variations 
were predicted as non-synonymous. In the Perm-R 
strain, ten genes were specifically affected, including 
two short-chain dehydrogenases, two heme peroxidases, 
two prophenoloxydases, one alcohol dehydrogenase, the 
glutathione S-transferase GSTE4 and one ABC trans- 
porter. Ten variations affecting seven genes including 
the GSTDl, three glycosyltransferases, one aldehyde 
oxidase, one ABC transporter and one oxidase/peroxid- 
ase were specific to the Imida-R strain. Finally, 24 varia- 
tions affecting seven genes including the GSTDl, one 
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Figure 5 Focus on detoxification: cytochrome P450s. Clustering 
analyses of transcription level variations (upper panel) and 
differential SNPs (lower panel) affecting cytochrome P450 
monooxygenases (CYP genes). Clustering was performed on all 
transcripts showing significant transcription level variations or 
differential SNPs within their coding region in any strain. The 
green-red color scale indicates transcription level variations as 
compared to the susceptible strain. Stars indicate a significant 
differential transcription. Blue-yellow color scale indicates allele 
frequency variations as compared to the susceptible strain. Amino 
acid position, amino-acid change, transcript number and gene 
names are indicated. Non-synonymous variations according to 
genome annotation are underlined. SNPs falling within P450 
Substrate Recognition Sites (SRS) are indicated. Names ending with 
a question mark indicate genes with ambiguous gene name 
(subfamily indicated). 



glycosyltransferase, one heme peroxidase, one carboxy- 
lesterase, two ABC transporters and one oxidase/perox- 
idase were specific to the Propo-R strain. 

Discussion 

Characterizing resistance mechanisms is essential for 
improving resistance management. Although target site 
modifications play a major role in resistance [4,5,28,29], 
other mechanisms such as insecticide biodegradation, al- 
tered transport, sequestration and modification of the 
insect cuticle also account for a significant part of resist- 
ance [6,18,30,31]. However, the intricacy of these mecha- 
nisms makes challenging the identification of candidate 
genes for functional validation [21]. Gene expression mi- 
croarrays are mostly used for identifying genes differen- 
tially transcribed in resistant populations but suffer from 
technical biases [21,32,33]. In contrast RNA-seq generates 
transcription data with a higher resolution, better dynamic 
range and lower technical variation [27]. In addition, 
RNA-seq produces polymorphism data and useful infor- 
mation to re-annotate gene models [25,26]. Despite recent 
studies pointing out the role of polymorphism variations 
in insecticide resistance [22-24], such aspect has never 
been investigated at the transcriptome level in mosquitoes. 
Indeed, the only study using RNA-seq to investigate in- 
secticide resistance in mosquitoes did not consider poly- 
morphism variations [34] . In this regard, the present study 
represents the first attempt to use RNA-seq for examining 
concomitantiy quantitative and qualitative transcriptome 
changes associated with resistance to different insecticides 
in mosquitoes. 

Insecticide selection and resistance levels 

After 10 generations of selection, bioassays revealed a 
constitutive increased resistance of each selected strain 
to its respective insecticide compared to the susceptible 
parental strain. Although resistance levels were low as 
compared to what can be observed in natura, they were 
significant regarding the few generations of selection. 
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Perm-R Imida-R Propo-R 




Perm-R Imida-R Propo-R 




■RB aldo-keto reductase 
■RA aldehyde dehydrogenase 
■RB glualhione S-transferase GSTD1 
-RA glucosyl /glucoronyl transferase 
-RB glucosyl /glucoronyl transferase 
-RA short-chain dehydrogenase 
-RA glucosyl / glucoronyl transferase 
-RA ABC transporter 
-RA glutathione S-transferase GSTD1 
■RC glutathione S-transferase GSTD1 
■RA glutathione S-transferase GSTE4 
i-RA cytochrome b5 
-RA glutathione S-transferase GSTD4 
-RA glucosyl /glucoronyl transferase 
-RA aldehyde oxidase 
-RA alcohol dehydrogenase 
■RA ABC transporter 
-RA glucosyl /glucoronyl transferase 
■RA glucosyl /glucoronyl transferase 
■RA heme peroxidase HPX5 
■RA oxidase/peroxidase 
-RA multicopper oxidase 
-RA heme peroxidase HPX1 
-RA prop henolox Ida se PP06 
-RA oxidase/peroxidase 
■RA prophenoloxidase PP09 
■RA Ihioredoxin peroxidase TPX3 
■RA short-chain dehydrogenase 
-RA lactoglutathione lyase 
■RD aldo-keto reductase 
■RC aldo-keto reductase 
-RA aldo-keto reducatse 
-RA short-chain dehydrogenase 
-RA alcohol dehydrogenase 
-RA s u If otransf erase 
■RA glucosyl / glucoronyl transferase 
■RA prophenoloxidase PP03 
-RA disulfite oxidoreductase 
-RA aldehyde oxidase 
-RA ABC transporter 
■RA su If otransf erase 
■RA heme peroxidase HPX2 
-RB ABC transporter 



261 P>P AAEL002354-RA heme peroxidase HPX5 
70 S>S AAEL001 586-RA glucosyl /glucoronyl tra 
70 S>S AAEL001586-RBglucosyl /glucoronyl transfei 
160 l>l AAEL001586-RA glucosyl /glucoronyl transfer 
160 l>l AAEL001586-RB glucosyl /glucoronyl trs " 
102 M>L AAEL001586-RA glucosyl /glucoronyl tr 
102 M>L AAEL001586-RB glucosyl /glucoronyl tr 
100 R>R AAEL001586-RAglucosyl-glucoronyltransferase 
100 R>R AAEL001586-RB glucosyl /glucoronyl transferase 
3 F>F AAEL001 061 -RA glutathione S-transferase GSTD1 
3 F>F /\AEL001061-RB glutathione S-transferase GSTD1 
3 F>F AAEL001061-RC glutathione S-transferase GSTD1 
1 1 P>P AAEL001 061-RA glutathione S-transferase GSTD1 
11 P>P AAEL001061-RB glutathione S-transferase GSTD1 
11 P>P AAEL001061-RC glutathione S-transferase GSTD1 
120 Q>Q /\AEL001061-RC glutathione S-transferase GSTD1 
125 AAELOOiOei-RCgiutathioneS-tfansferaseGSTDI 

189 E>E AAELOOIOei-RBglutathioneS-transferaseGSTDI 

190 A>A AAELOOIOei-RBglutathione S-transferase GSTD1 
965 R>R /\AEL001804-FiA glucosyl /glucoronyl transferase 
1513 G>G AAEL012698-RA ABC transporter 

298 V>V AAEL011763-RAprophenolcxidasePPO3 
729 L>L AAEL012698-RA ABC transporter 
806 >\ AAEL012698-FiA ABC transporter 
428 L>L AAEL014837-F!A prophenoloxidase PP09 
263 l>i AAEL0U371-RA glucosyl /glucoronyl transferase 
1 387 G>G AAEL012698-RA ABC transporter 
388 E>E AAEL009948-RA aldehyde dehydrogenasae 
329 S>S AAEL013171-FiA heme peroxidase HPX2 
51 G>G /\AEL014935-RA cytochrome bS 
352 y>L AAEL002016-RAdisultite oxidoreductase 
278 E>E AAEL00965-RA short-chain dehydrogenase 
280 A>A AAEL009625-RA short-chain dehydrogenase 
326 R>R AAEL002683-RA aldehyde oxidase 
199 R>R /\AEL005416-RA oxidase/peroxidase 
425 F>F AAEL002683-FiA aldehyde oxidase 
119 F>F AAEL014548-RAthiore(loxin peroxidase TPX3 
127 T>T AAEL014548-F!Athioredcxin peroxidase TPX3 
470 U>D AAEL012192-RA ABC transporter 
225 P>P AAEL014244-RA glucosyl /glucoronyl transferase 
1071 l>l AAEL012698-RAABCtransnortor 
zj4 T>T AAEL002354-RA heme peroxidase HPX5 
.87 A>A AAEL002354-RA heme peroxidase HPX5 
450 L>L AAEL0D2354-RA heme peroxidase HPX5 
453 K>K AAEL002354-FiA heme peroxidase HPX5 
1085 D>D /\AEL004968-FiA ABC transporter 
248 P>P AAEL002354-F!A heme peroxidase HPX5 
1 1 L>L AAEL007952-RA glutathione S-transferase GSTE4 
625 D>A AAEL005416-RAoxidase/peroxidase 
268 L>L AAEL005416-RA oxidase peroxidase 
1201 E>E AAEL012698-FiA ABC transporter 
272 A>S AAEL003195-RAcarboxi/choline esterase CCEAE1C 
619 A>A AAEL005416-RA oxidase/peroxidase 
292 H>H AAEL005416-RA oxidase/peroxidase 
28 l>l AAEL013171-RA heme peroxidase HPX2 
122 E>I AAEL010386-RA glucosyl /glucoronyl transferase 
339 G>G AAEL005416-F!A oxidase/peroxidase 
414 R>R AAEL014244-RA glucosyl /glucoronyl transferase 
183 MH AAEL002016-RAdisulfiteoxidoreductase 
107 K>K AAEL011130-RAalcohol dehydrogenase 
112 L>L AAEL011130-RAalcoholdehydrogenase 
536 R>R AAEL002016-RAdisulfteoxidoreducatse 
140 L>L AAEL011130-RAalcoholdehydrogenase 
651 Q>Q AAEL002016-RAdisulfite oxidoreductase 
199 R>R AAEL012401-RAshort-chian dehydrogenase 
284 A>y AAEL009625-RA short-chain dehydrogenase 



Figure 6 Focus on detoxification: Other enzymes and 
transporters. Clustering analyses of transcription level variations 
(upper panel) and SNPs (lov\/er panel) affecting enzymes and 
transporters potentially involved in Insecticide detoxification. Only 
transcripts or differential SNPs falling within coding regions are 
shown. Green-red color scale indicates transcription level variations 
as compared to the susceptible strain. Stars indicate a significant 
differential transcription. Blue-yellow color scale indicates allele 
frequency difference as compared to the susceptible strain. Amino 
acid position, amino-acid change, transcript number and gene 
names are indicated. Non-synonymous variations according to 
genome annotation are underlined. 



the full susceptibility of the parental strain and the ab- 
sence of target-site resistance alleles. The rapid rise of 
resistance suggests that alleles conferring a better fitness 
in presence of insecticides are already present in suscep- 
tible populations and can be promptly selected under a 
constant selection pressure. Although laboratory selection 
does not fully mimic adaptive processes occurring in nat- 
ura (e.g. lower population sizes, different environment, no 
introduction of resistant alleles by migration ...), selecting 
all resistant strains from a single susceptible strain allowed 
minimizing variations related to different genetic back- 
grounds. Using a fully susceptible parental strain also 
allowed focusing on non-target site resistance mecha- 
nisms (no target-site mutations detected in parental or se- 
lected strains). 

RNA-seq as a tool for studying insecticide resistance in 
mosquitoes 

A total of 270 million cDNA reads were sequenced across 
all samples and 80% of them were successfully mapped to 
Ae. aegypti reference genome. Such mapping efficiency ap- 
peared acceptable considering polymorphism variations 
occurring between the reference genome (Liverpool strain) 
and the parental strain used in our study (Bora-Bora 
strain). The high correlation of expression data obtained 
from cDNA library replicates confirmed the robustness of 
cDNA library construction and sequencing procedures. By 
applying high-stringency sequence quality and transcript 
coverage filtering, high fidelity transcription data were re- 
covered for more than 13000 transcripts. Such detection 
level was comparable to those obtained with DNA micro- 
arrays at the same life stage [35,36] . More than 2500 differ- 
ential SNPs linked to insecticide selection were identified. 
Although our experimental design did not control for sto- 
chastic effects (genetic drift) and the presence of false posi- 
tives is likely, transcripts affected by these differential SNPs 
represent strong candidates for further fiinctional valid- 
ation studies. Little overlap was found between transcripts 
differentially transcribed in selected strains and those af- 
fected by differential SNPs. This was expected as RNA-seq 
data are restricted to transcripts and did not cover regula- 
tory regions often located outside transcript boundaries. 



David ef al. BMC Genomics 2014, 15:174 
http://www.biomedcentral.com/1471-2164/15/174 



Page 9 of 15 



Challenging reads distribution with genome annotation 
identified > 3000 transcripts incorrectly annotated with 
most of them showing wrong UTR boundaries or modi- 
fication of exon/intron structure. In addition, more than 
500 lonely genomic regions showing high transcription 
signals and realistic exon/intron structures were identi- 
fied. Further analyses are now required for assigning 
them to new exons of known transcripts, novel tran- 
scripts, pseudogenes or non-coding RNAs. 

Transcriptome changes associated with insecticide 
selection 

Overall, our study revealed diverse response patterns de- 
pending on the insecticide used for selection. The Imida-R 
strain selected with the neonicotinoid imidacloprid showed 
the highest resistance level together with the strongest dif- 
ferential transcription response with numerous genes being 
over transcribed. In contrast, the Perm-R strain selected 
with the pyrethroid permethrin revealed a moderate differ- 
ential transcription response but a lot of polymorphism 
variations following insecticide selection. The Propo-R 
strain selected with the carbamate propoxur showed an 
intermediate pattern. These dissimilar quantitative and 
qualitative transcriptome changes may reflect different 
adaptive strategies driven by costs and benefits associ- 
ated with resistance mechanisms to each insecticide. 

Insecticide resistance and immunity 

Multiple transcripts related to immune response were 
affected by insecticide selection with several of them 
under transcribed in all resistant strains (cecropins, defen- 
sins, lectins, ...) and others affected by differential SNPs 
(defensins, clip-domain proteases, spatzle proteins...). 
Mosquito humoral response is involved in their capacity 
to host and transmit viruses [37] or parasites [38]. As 
mentioned in Alout et al. [39], our results support an as- 
sociation between insecticide resistance and the capacity 
of mosquitoes to host and transmit pathogens, which may 
affect the control of mosquito-borne diseases [40]. 

Altered insecticide transport and sequestration 

Multiple transcripts encoding cuticular proteins were 
over transcribed in the Imida-R strain while most of 
them were under transcribed in the Perm-R strain and 
not affected in the Propo-R strain. Such strong effect in 
the Imida-R strain was not associated to changes in the 
polymorphism of these transcripts, suggesting that insecti- 
cide selection is affecting their regulation rather than their 
conformation. Cuticle plays a crucial role in protecting in- 
sects from their environment. The vast majority of chem- 
ical insecticides are active by contact and changes in 
cuticle thickness or conformation have been suggested to 
contribute to resistance in mosquitoes [31,41,42]. Our 
data support previous studies suggesting a significant role 



of cuticle in the adaptation to neonicotinoid insecticides 
[36,43,44]. The role of cuticle in response to imidacloprid 
selection was supported by the specific over transcription 
of the multi-copper oxidase AAEL007415 in the Imida-R 
strain as its An. gambiae orthologue is involved in cuticle 
and egg shell tanning [45,46]. Numerous protein families 
are involved in cuticle biosynthesis and homeostasis in- 
cluding enzymes, transporters and transcription factors 
and further studies are now required for pinpointing those 
controlling cuticular resistance in mosquitoes. 

Several kinases were specifically over-transcribed in the 
Imida-R strain. Kinases are involved in multiple regulatory 
mechanisms and their involvement in the response of 
insects to insecticides is likely. Indeed, recent studies 
showed that the phosphorylation state of acetylcholine 
nicotinic receptors can modulate the efficacy of neonicoti- 
noid insecticides [47-49] . 

Five transcripts encoding hexamerins were strongly 
over transcribed in the Perm-R strain and in a lesser ex- 
tend in other resistant strains. These transcripts, located 
on different supercontigs, were also affected by differen- 
tial SNPs suggesting a selection imprint on these pro- 
teins. Insect hexamerins may be involved in hormone 
transport, energy and amino acid storage, cuticle biosyn- 
thesis and immune defense [50,51]. Hexamerins of the 
lepidopteran Heliothis zea have been shown to bind insec- 
ticides, suggesting a direct role in resistance by sequestra- 
tion or altered transport [52]. Deciphering if hexamerins 
are impacted by insecticide selection because of their abil- 
ity to bind insecticides, their interaction with cuticle 
homeostasis or their role in resources re-allocation associ- 
ated to fitness costs remains unclear. 

Finally, ATP-binding cassette transporters (ABC trans- 
porters) can play a role in adaptation to xenobiotics [53] 
and have been associated to insecticide resistance in mos- 
quitoes [30,54]. In our study, the response of ABC trans- 
porters to insecticide selection was marginal. Indeed, only 
one ABC transporter (AAEL008624) was found differen- 
tially transcribed in response to insecticide selection and 
this gene was down regulated in all resistant strains. Four 
others were affected by differential SNPs but their allele 
frequency variations in resistance strains were low. 

Insecticide biodegradation 

Detoxification enzymes play a major role in insecticide 
resistance [4,6,13,14,18,55]. As expected, these enzymes 
were well represented in our data set but showed dis- 
tinct patterns depending on the nature of the insecticide 
used for selection. 

Response to selection with the neonicotinoid imidaclo- 
prid (Imida-R strain) was characterized by the over tran- 
scription of multiple P450s, oxidases, transferases and one 
alcohol dehydrogenase, supporting the involvement of mul- 
tiple enzymes in imidacloprid biodegradation pathways. 
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Among them, the P450 CYP6BB2 was recently pointed out 
as a solid candidate for imidacloprid metabolism based on 
gene expression data and substrate binding predictions 
[36]. The aldehyde oxidase AAEL002683 was over tran- 
scribed and affected by polymorphism variations in the 
Imida-R strain. A recent study confirmed that aldehyde oxi- 
dases can contribute to neonicotinoid metabolism through 
nitro-reduction [56]. 

Several P450s over transcribed in the Imida-R strain 
were also found over transcribed in response to propo- 
xur selection. Cross resistance between these two insec- 
ticides was identified [36] and the potency of particular 
P450s to confer resistance to multiple insecticides has 
been shown. For example, An. gambiae CYP6Z1 metab- 
olizes the organochlorine DDT and the carbamate car- 
baryl [22] while CYP6M2 metabolizes both DDT and 
pyrethroids [57,58]. 

Metabolic resistance of mosquitoes to pyrethroids 
has been mainly associated with an over expression of 
P450s able to metabolize them. Among the multiple 
candidates identified by microarray screenings, An. gam- 
biae CYP6M2 and CYP6P3, An. funestus CYP6P9, Ae. 
aegypti CYP9J32, CYP9J24, CYP9J28, and Cx quinque- 
fasciatus CYP6M10 and CYP4H24 have been validated 
as pyrethroid metabolizers [review in 18]. CYP6Zs and 
CYP6Ms have also been associated with pyrethroid re- 
sistance by QTL in An. funestus where resistance mainly 
relies on metabolic mechanisms [59]. Recentiy, the central 
role of mosquito CYP6Zs in pyrethroid degradation path- 
way was revealed [55]. In the present study few P450s were 
found over transcribed in the Perm-R strain. Among them, 
CYP6BB2 was also found strongly over transcribed in pyr- 
ethroid resistant populations from Cuba and Cayman 
islands [30]. Unexpectedly, known Ae. aegypti pyrethroid 
metabolizers or paralogs of those validated from other 
mosquito species were not found over transcribed in the 
Perm-R strain. However, a strong selection imprint was de- 
tected in the Perm-R strain in a P450 cluster in supercontig 
1.371containing strong candidates CYP6Ms, CYP6Ns and 
CYP6Zs, suggesting that the selection of particular P450 
variants can contribute to pyrethroid resistance. In mam- 
mals, it is well known that P450 variants can display differ- 
ent substrate specificity [60-64]. To date, such variations 
have been neglected in mosquitoes with only few studies 
pointing out P450 alleles associated with insecticide resist- 
ance [22,23]. Finally, the strong selection imprint observed 
in these P450s might explain their unexpected under tran- 
scription in the Perm-R strain. Indeed, the apparent under 
expression of transcripts affected by differential SNPs 
might be the consequence of a mapping bias due to a 
higher divergence from the reference genome or an enrich- 
ment in low-expressed alleles [65,66]. Although further 
analyses are required for investigating the role of allele- 
specific expression in resistance, our data supports the 



selection of particular detoxification enzyme variants by 
insecticides. 

Conclusions 

The present study primarily aimed at assessing the use- 
fulness of RNA-seq to investigate insecticide resistance 
mechanisms in mosquitoes. Results confirmed that this 
technique produces high-quality gene expression data to- 
gether with solid polymorphism data. Distinct responses 
to selection with insecticides from different chemical 
families were observed with a balance between gene 
expression and polymorphism variations. Polymorphism 
variations of P450 enzymes were strongly linked to pyreth- 
roid selection. Although additional analyses are required 
to validate variants linked to resistance, such finding high- 
lights the necessity to consider both gene expression and 
polymorphism variations for identifying candidate genes 
potentially involved in insecticide resistance. As sequen- 
cing costs are decreasing and new sequencing strategies 
and bioinformatics pipelines are developed, obtaining gene 
expression and polymorphism data from the same sam- 
ples using high-throughput sequencing should now be 
considered as a valuable alternative to microarrays. 

Methods 

Mosquito selection with insecticides and bioassays 

The mosquito Ae. aegypti, was used in the present study. 
Mosquitoes were reared in standard insectary conditions 
(26°C, 14 h/10 h light/dark, 80% relative humidity) in 
tap water (larvae) and net cages (adults). Larvae and 
adults were fed with hay pellets and papers impregnated 
with honey respectively. Blood feeding of adult females 
was performed on mice. Mice were maintained in an 
animal house agreed by French Ministry of animal wel- 
fare (n° B 38 421 10 001) and used in accordance to EU 
laws and the relevant ethic committee recommendations 
(ComEth Grenoble - C2EA - 12). The laboratory strain 
Bora-Bora, originating from French Polynesia and fully 
susceptible to insecticides, was used as a parental strain 
to select three independent strains with the pyrethroid 
insecticide permethrin (Perm-R strain), the neonicoti- 
noid insecticide imidacloprid (Imida-R strain) and the 
carbamate insecticide propoxur (Propo-R strain). Both 
pyrethroid and carbamate insecticides are heavily used 
against mosquitoes and resistance to these insecticides 
has been reported worldwide. Neonicotinoids are mar- 
ginaly used against mosquitoes but represent one possible 
alternative when resistance to other insecticide threatens 
the efficacy of mosquito control [36]. Selection was per- 
formed by exposing early 4th-stage larvae for 24 h to a le- 
thal dose of each insecticide. For each strain, the dose of 
insecticide was adjusted at each generation (4 to 6 |ig/L 
permethrin, 500 to 1000 [ig/L imidacloprid and 500 to 
800 ng/L propoxur) in order to reach 60-80% larval 
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mortality after 24 h exposure. Surviving larvae were trans- 
ferred in clean tap water, fed with standard larval food and 
allowed to emerge. Adults were allowed to reproduce for 
4 days and blood-fed to obtain eggs for the next gener- 
ation. In order to limit bottleneck effects, each generation 
was seeded with more than 6000 individuals. Selection 
process was carried out in parallel for all strains during 
ten generations. During this process, the susceptible par- 
ental strain was maintained in similar condition without 
insecticide selection. Bioassays and molecular work were 
performed on early 4th-stage larvae of the 11th generation 
(Gil larvae) bred in standard conditions and not exposed 
to insecticides. 

To assess the constitutive resistance level of each se- 
lected strain, larval bioassays were conducted with per- 
methrin, imidacloprid and propoxur comparatively to the 
susceptible parental strain. Five doses of each insecticide 
and four replicates of 25 larvae per dose were used. Doses 
of permethrin (1.5 to 6.5 \xg/V), imidacloprid (150 to 
2200 i^g/L) and propoxur (100 to 2000 |ig/L) were chosen 
in order to cover the whole mortality range after 24 h 
exposure. Lethal concentrations corresponding to 50% 
mortality (LC50) and their 95% confident intervals (CI950/J 
were then calculated with a probit approach for each 
strain using XL-Stat (Addinsoft, Paris, France). Resistance 
ratios (RR50 based on LC50 values) were calculated by 
comparison to the susceptible parental strain. 

RNA extraction and cDNA libraries preparation 

For each strain, total RNA was extracted from Gn larvae 
obtained from three independent egg batches (three 
biological replicates per strain). Each biological replicate 
consisted of 180 larvae reared in 200 mL tap water in 
standardized insectary conditions. For each biological 
replicate, total RNA was extracted from 60 4th-stage lar- 
vae using the RNAqueous-4PCR kit (Applied Biosystems/ 
Ambion, Austin, TX, USA), according to the manufac- 
turer's instructions. Total RNA quality and quantity were 
assessed with a Nanodrop NDIOOO (Thermo Scientific, 
USA) and a 2100 Bioanalyzer (Agilent, USA). After extrac- 
tion, total RNAs from each biological replicate were 
pooled in equal quantity to obtain a total RNA mixture 
representative of 180 individuals. Total RNA pools from 
each strain were then used for preparing cDNA libraries 
using the mRNA-Seq Sample Prep Kit (Part 1004898 Rev 
D, lUumina, USA). Two replicates of cDNA libraries were 
prepared for each strain as follows. Briefly, mRNAs were 
purified using poly-T beads and chemically fragmented. 
These fragmented mRNAs were reverse-transcribed using 
Superscript II (Invitrogen) at 42°C for 50 min. Double- 
stranded cDNAs were then synthesized and mRNAs re- 
moved using DNA polymerase I and RNase H at 16°C for 
2.5 hours. Double-stranded cDNAs were purified using 
QIAquick PGR Purification Kit (QIAGEN, Germany) and 



processed for end-repair and 3 ' adenylation using Klenow 
polymerase. Sequencing adaptors were then ligated using 
DNA ligase. Adapter-ligated cDNA libraries were then 
purified on 2% agarose gel based on a size range of 200 ± 
25 bp. Adapter-ligated cDNA libraries were then enriched 
by 15 PGR cycles with adaptor-specific primers using 
Phusion DNA polymerase (Finnzymes Oy). Enriched 
cDNA libraries were then purified using QIAquick PGR 
Purification Kit (QIAGEN, Germany) before quality 
control analysis on an Agilent 2100 Bioanalyzer. 

Sequencing, read mapping and genome re-annotation 

Each double-stranded cDNA library was sequenced as 
single reads of 75 bp in a distinct flow cell lane with a 
Genome Analyzer II (lUumina) at the National Sequen- 
cing Center (Genoscope, Evry, France). Reads were then 
mapped to the Ae.aegypti genome sequence (Aaeg LI. 2 
gene set, http://vectorbase.org) using the Tophat algo- 
rithm with default parameters (http://tophat.cbcb.umd. 
edu, release 1.0.14) [67], leading to an average mapping 
rate of 80%. Bam files were then loaded into Genespring 
NGS version 12.5 (Agilent Technologies) for further ana- 
lysis. Reads were then filtered based on the sequence qual- 
ity (mean read quality > 30 and with < 10 Ns) and mapping 
quality (alignment score > 98), and non-primary multiply 
mapped reads were removed. The remaining reads were 
used for all subsequent analyses. Read coverage was com- 
puted and genome annotation was challenged based on 
read distribution and coverage. Default feature detection 
parameters were used (min exon length percentile 10; min 
intro length percentile 10; max intron length percentile 90; 
min exon RPKM percentile 50; min gene RPKM percentile 
50; min gene length percentile 10; min exon RPKM with 
respect to host gene percentage 75; min number of reads 
in exon 10). Modified exon-intron structures and new pu- 
tative transcripts were proposed based on read coverage, 
read splicing events and distance with respect to existing 
genes and transcripts. Reads per KUobase exon Model per 
Million sequenced reads (RPKMs) were calculated for all 
known and putative new transcripts. Transcript RPKM 
values were used for assessing the variability between li- 
brary replicates for each strain. 

Differential transcription analysis 

In order to avoid estimating transcription ratios from low 
coverage transcripts, only the 13105 transcripts showing 
at least 0.5 RPKM per condition (~ 30 reads/kb) were con- 
sidered for differential transcription analysis. Differential 
transcription of each transcript was then tested between 
insecticide-selected strains and the susceptible parental 
strain using an Audic-Claverie test (AG test) based on read 
counts [68] with a Benjamini and Hochberg's multiple 
testing correction [69]. Transcripts showing an adjusted 
P value < 10"^^ and a fold change > 3 in either direction 
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were considered as differentially transcribed between 
insecticide-selected strains and the parental susceptible 
strain. RNA-seq transcription data from the Imida-R 
strain were compared with transcription data obtained 
from the same biological samples with the Aedes aegypti 
15 K DNA microarray as previously described [36]. Only 
the 326 transcripts showing significant adjusted P values 
in both techniques were considered. 

Polymorphism detection and differential analysis 

Detection of polymorphisms was performed based on 
the 174 million reads passing quality filters with the fol- 
lowing parameters: confidence score threshold = 100, 
coverage > 20 reads, base quality cut off = 5, ignore loca- 
tions within or next to homopolymer stretches > 10 nu- 
cleotides. Among all detected polymorphism variations, 
only SNP substitutions were considered for differential 
polymorphism analyses. SNP allele frequencies were then 
computed between each insecticide-selected strain and 
the susceptible strain. Allele frequencies were considered 
as differential between an insecticide-selected strain and 
the susceptible strain (hereafter named as differential 
SNPs) if the following conditions were fulfilled: i) Total 
read coverage at SNP position between both strains > 50, 
ii) Strand bias at SNP position < 50 for both strains and 
Hi) Allele frequency difference between both strains > 50% 
in either direction. Potential genie effects of SNPs were 
computed by comparing SNP with reference genome 
annotation. Genie effects were defined as 5'UTR, Syn- 
onymous, Non-synonymous, Intronic, 3 'UTR and Inter- 
genic {i.e. close but not within gene boundaries). 

GO term enrichment analyses 

Go term enrichment analyses were conducted on i) tran- 
scripts significantly over transcribed in insecticide-selected 
strains compared to the susceptible strain, ii) transcripts 
significantly under transcribed in insecticide-selected 
strains compared to the susceptible strain and Hi) tran- 
scripts affected by differential SNPs in their coding se- 
quence as compared to the susceptible strain. GO term 
enrichment analyses were performed separately on each 
transcript list versus all detected transcripts. Enrichment 
statistics were computed based on hypergeometric distri- 
bution and adjusted with Benjamini-Yekutieli's multiple 
testing correction [70] which takes into account the 
dependency among GO terms. GO terms showing an ad- 
justed P value < 0.05 were considered as significantly 
enriched in insecticide-selected strains compared to the 
susceptible strain. 

Clustering analyses 

All clustering analyses were performed with TM4 MEV 
version 4.3.02 [71] using Euclidean distance and complete 
linkage algorithm with optimization of gene and condition 



trees. Transcripts significantly differentially transcribed 
(see above for thresholds) in any insecticide-selected strain 
compared to the parental susceptible strain were clustered 
based on their TR (log2 ratios). Annotated transcripts 
represented in the main clusters were then assigned to 
nine biological categories based on their annotation (in- 
cluding GO terms). Categories were defined as follow: 
'Detoxification', 'Kinases/Phosphatases', 'Other enzymes', 
'Cuticle', 'Immunity', 'Hormones and neurotransmitters 
signaling', 'Transcription factors', 'Intra or extracellular 
trafficking/chaperonins' and 'Structure'. Enrichment of 
these categories compared to all detected transcripts was 
then computed using a one-side Fisher's exact test 
followed by Benjamini and Hochberg multiple testing cor- 
rection [69]. Categories showing a corrected P value < 0.05 
were considered significandy enriched. A similar cluster- 
ing analysis was performed with differential SNPs (see 
above for detailed criteria) falling in coding regions. Clus- 
tering was performed as described above using allele fre- 
quency variations (-100% to +100%) in each strain as 
compared to the susceptible strain. Differential SNPs were 
assigned to nine different biological categories (see above) 
in regard of the nature of the transcript affected. For each 
cluster, biological category enrichment was computed ver- 
sus all detected transcripts as described above. 

Clustering analyses were then focused on transcripts en- 
coding genes potentially involved in insecticide detoxifica- 
tion pathways (detoxification sensus lata). Clustering of 
both transcription level variations and differential SNPs 
were performed as described above. Separate analyses 
were conducted for cytochrome P450 monoooxygenases 
(CYPs) and other detoxification transcripts. 

Additional files 

r ^ 

Additional file 1: Table SI. Sequencing and mapping statistics. 

Additional file 2: Figure SI. RPKM correlation between cDNA library 
replicates. Each dot represents one transcript Only transcripts showing 
more than 0.5 RPKM are shown. 

Additional file 3: Figure S2. Comparison of read coverage across 
strains. Read coverage are indicated for each strain as RPKM (log scale). 
Coverage distributions are compared for unmodified transcripts (top), 
re-annotated transcripts (bottom left), and new putative transcripts 
(bottom right). 

Additional file 4: Table S2. Overview of transcriptome re-annotation. 

Additional file 5: Table S3. Genomic location of all novel transcribed 
features. 

Additional file 6: Figure S3. Cross-validation of transcription levels 
between RNA-seq and microarrays. Comparison is based on transcription 
data obtained from the Imida-R strain versus susceptible strain. RNA-seq 
and microarray data were obtained from the same generation. 
Correlation was performed on the 326 transcripts showing a significant 
differential transcription level in both studies. Solid grey line represents 
an equal transcription ratio between both techniques. Grey dashed lines 
represent a two-fold variation. 

Additional file 7: Table S4. Transcripts differentially expressed after 
insecticide selection. For each transcript transcription level (fold) 
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compared to the susceptible strain, adjusted P value and RPKIVl are 
indicated. Clusters and biological functions as described in Figure 2 are 
indicated. 

Additional file 8: Table S5. Differential SNPs linked to insecticide 
selection. For each SNP, allele frequency variations in each strain 
compared to the susceptible strain are indicated together with the 
affected transcript the cDNA position, the amino-acid position and the 
amino acid change as compared to the reference genome. Clusters and 
biological functions as described in Figure 4 are indicated. 

Additional file 9: Figure S4. GO terms enrichment analysis. Analysis 
was performed on all transcripts significantly differentially expressed or 
affected by differential SNPs in insecticide-selected strains as compared 
to the susceptible strain. GO terms associated to each transcript were 
extracted from Vectorbase. GO terms showing adjusted P values < 0.05 
were considered significantly enriched. 

Additional file 10: Figure S5. Differential SNPs linked to permethrin 
selection in supercontig 1.371. Transcripts location and read coverage are 
indicated. Transcripts showing differential SNPs in the Perm-R strain as 
compared to the susceptible strain are shown in red. For each transcript, 
the number of differential SNPs and their predicted genie effects are 
indicated. 
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